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ABSTRACT 

We determine the mass of the black hole at the center of the spiral galaxy NGC 
4258 by constructing axisymmetric dynamical models of the galaxy. These mod- 
els are constrained by high spatial resolution imaging and long-slit spectroscopy 
of the nuclear region obtained with the Hubble Space Telescope, complemented 
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by ground-based observations extending to larger radii. Our best mass estimate 
is M. = (3.3 ± 0.2) x 1O 7 M for a distance of 7.28 Mpc (statistical errors only). 
This is within 15% of (3.82 ± 0.01) x 1O 7 M , the mass determined from the 
kinematics of water masers (rescaled to the same distance) assuming they are 
in Keplerian rotation in a warped disk. The construction of accurate dynamical 
models of NGC 4258 is somewhat compromised by an unresolved active nucleus 
and color gradients, the latter caused by variations in the stellar population 
and/or obscuring dust. Depending on how these effects are treated, as well as on 
assumptions about the ellipticity and inclination of the galaxy, we obtain black 
hole masses ranging from 2.4 x 1O 7 M to 3.6 x 1O 7 M . This spread is mainly due 
to uncertainties in the stellar mass profile inside the central 2" (~70 pc). Ob- 
scuration of high- velocity stars by circumnuclear dust (possibly associated with 
the masing disk) could lead to an underestimate of the black hole mass which 
is hard to correct. These problems are not present in the ~ 30 other black hole 
mass determinations from stellar dynamics that have been published by us and 
other groups; thus, the relatively close agreement between the stellar dynamical 
mass and the maser mass in NGC 4258 enhances our confidence in the black hole 
masses determined in other galaxies from stellar dynamics using similar methods 
and data of comparable quality. 

Subject headings: galaxies: spiral — galaxies: individual (NGC 4258) — galaxies: 
kinematics and dynamics — galaxies: nuclei 

1. INTRODUCTION 

The characterization of massive dark objects (hereafter "black holes") at the centers of 
nearby galaxies has advanced considerably in the past decade. Nonetheless the dynamical 
modeling of both stellar and gas kinematics remains subject to uncertainties. These can 
arise from, among other things, systematic uncertainties and noise amplification in the de- 
projection of the observables from two to three dimensions, the presence of dust that can 
affect both the mass-to-light ratio (M/L) and the observed kinematics, variations in M/L 
due to stellar population gradients, and the presence of nonstellar point and/or continuum 
sources in the central galactic regions which may swamp the light from the surrounding stars 
or gas. Gas kinematics can also be affected by nongravitational forces, random velocities 
and pressure support, which are difficult to model accurately, while the complex phase-space 
structure of stellar orbits can make it difficult to properly sample initial conditions and the 
phase-space density distribution of the orbits in (stellar) galactic equilibria. Reverberation 
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mapping presents a third po ssibility for determin i ng black hole masses but comes with its 



own set of limitations (e.g., iBlandford fc McKed Il982l ; iPetersonl Il993l ; iNetzer fc Peterson 



19971 : iGebhardt et alJl2000bh 



It is therefore desirable to test the reliability of these methods. Unfortunately, com- 
paring stellar dynamical estimates to gas dynamical estimates of black hole masses is not 
usually possible, as few galaxies have properties that permit more than one method to be 
applied with confidence. Even in those cases where it is possible, a detailed comparison 
is not always illuminating due to uncertainties in the measurements. One case where gas 
dynamical and stellar dynamical mass measurements are available for the same galaxy is not 
reassuring: the gas dynamical m ass for the black hole in IC 1459 is a factor of ~ 5 smaller 
than the stellar dynamical mass (IVerdoes Kleijn et al.ll2000l : ICappellari et al.ll2002l ). Better 
agreement, within a factor of ~ 2, was found between the mass estimate of the black hole 



in Centaurus A by ISilge et al.l (120051 ) from stellar dynamics, and by iMarconi et al.l (120061 ) 
from gas dynamics. Even in the case of the supermassive black hole in t he nucleus of our 
own Milky Way Galaxy, the mass estimates from individual stellar orb its (jGhez et ah 2005 ) 
and f rom a statistical analysis of radial velocities and proper motions (jChakrabarty fc Saha 
20011 ) differ by about a factor of two. 



NGC 4258 presents a uniq ue op portunity to test the re liability of black hole mass es- 



timators. iMiyoshi et al.l (119951 ) and iHerrnstein et al.l (119991 ) used the Very Long Baseline 
Array (VLBA) to observe water maser emission in a thin, warped, and nearly edge-on gas 
annulus located at 0.16-0.28 pc from the center of this galaxy, and determined that a central 
black hole with a mass M. = (3.9 ±0.1) x 10 7 M Q is needed to account for the observed 
velocity profile, proper motions, and accelerations. This black hole mass is thought to be 
quite reliable for several reasons, (i) The maser emission lines are both strong and narrow, 
lending themselves to VLBA observations of much higher resolution than possible with the 
Hubble Space Telescope (HST) or ground-based telescopes, (ii) The masers seem to lie in 
a disk that is very thin (height /radius ~ 0.02), so the corrections to the rotation curve due 
to pressure or random velocities are negligible, (iii) Most importantly, the maser velocities 
scale very nearly as v (r) oc r -1 / 2 with distance r from the center, so that they exhibit an 
almost perfect Keplerian motion, which makes the deduction of the black hole mass easy 
and almost free of assumptions, (iv) Finally, the measurements of the mass and distance 
from the radial velocities, proper motio ns, and accelerations of the masers are all consistent. 
More recently, IHerrnstein et al.l (120051 ). using a maximum likelihood analysis of the maser 
positions and velocities, quantified the deviation from Keplerian motion. A number of dif- 
ferent models that can explain this discrepancy yield central masses in the range (3.59-3.88) 
xl0 7 M Q rescaled to a distance of 7.28 Mpc. The warped disk model, which the authors 
consider the most probable, implies M. = (3.82 ± 0.01) x 1O 7 M at this distance. 
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In this paper, we estimate the mass of the central black hole of NGC 4258 via mod- 
eling of the stellar dynamics of the galaxy. This is accomplished using the orbit super- 
position tools developed by our collabo r ation (the "Nuker" team) (IGebhardt et al.ll2000al ; 



Bower et all bpOll : bebhardt etaD l20oij~TIiomas et al 



tion tools have bee n developed by Ivan der Marel et al. 



2004T). Alternative orb it supe rposi- 



(Il998h . ICretton et all (11999 ). and 



Valluri et al.l (120041 ) . Because the maser black hole mass estimate seemed very accurate, we 
intended this to be a thorough, end-to-end test of the accuracy of our method, at least for 
this galaxy. Unfortunately, we encountered a number of problems that hindered our ability 
to determine unambiguously the mass density profile in the central region of the galaxy. In 
particular, (i) it was not possible to accurately subtract the light of the unresolved active 
galactic nucleus (AGN) at the center of the galaxy to recover the underlying stellar luminos- 
ity profile; (ii) it was not possible to reliably determine the stellar M/L profile within ~2" 
of the nucleus due to a color gradient, presumably caused by diffuse dust and/or a spatial 
variation in the stellar population; and (iii) it was not possible to estimate the extinction 
caused by a compact dust patch near the nucleus, which may be obscuring light from fast- 
moving stars very close to the black hole, thus causing a systematic underestimate of the 
black hole mass. 

Despite these difficulties, our "most trusted" estimate, M. = (3.3 ± 0.2) x 10 7 M Q , is in 
reasonable agreement with the maser mass estimate. The quoted error margin corresponds 
to the "la" confidence band, and reflects statistical errors only. Uncertainties in the mass 
profile discussed in the previous paragraph widen the error margin to 2.4 x 1O 7 M < M. < 
3.6 x 1O 7 M . This range still does not include the potential effect of obscuring circumnuclear 
dust, which cannot be easily estimated, or any other systematic effects such as deviations 
from axisymmetry that cannot be properly treated by our dynamical modeling method. 
These factors are discussed in detail in the following sections. 



Pastorini et al.l (120071 ) estimate the mass of the black hole in NGC 4258 using gas kine- 
matics, but their 95% confidence intervals span a factor of ten in mass so this is not a strong 
test of the reliability of black-hole mass estimates. 



For this work, we use 7.28 Mpc for the distance to NGC 4258 (ITonry et al.ll200ll ). which 
is in excellent agreement with the maser-derived distance of 7.2 Mpc. At that distance, 
1" ~ 35 pc. 

In £j2]we present the photometric and kinematic datasets, which consist of HST WFPC2 
and NICMOS imaging and STIS spectroscopy, complemented by ground-based imaging and 
long-slit spectroscopy that extend data coverage to larger radii. We discuss the large-scale 
morphology of the galaxy in §|3l mainly as it pertains to the determination of the contribution 
of the stars to the gravitational potential. In §H we consider the consequences of the color 
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gradient in the central 2". §5] presents the dynamical models. We conclude in §61 



2. OBSERVATIONS AND REDUCTION 
2.1. Photometry 

Here we describe the acquisition and reduction of the photometric data, which are 
needed mainly to determine the stellar mass distribution and the kinematic tracer profile, in 
the way described in §5.11 and §5.41 

A journal of the imaging observations, from which the surface photometry was derived, 
is shown in Table [TJ A wide variety of data are available in different colors, both from HST 
and from the ground. For reasons explained in §4.11 we produced three luminosity profiles, 
(i) The near-infrared (NIR) profile was created by combining ground-based and HST J, H, 
and .fT-band images, as described in §2.1.11 (ii) The V^-band profile was generated from 
images in the V band (HST) and in the R band (ground-based), as described in §2.1.21 
(iii) The J-band profile was created from the HST /-band image, as described in §2.1.21 
Although this image was saturated in the nucleus, it was useful to consult for the luminosity 
weighting (normalization) of the line-of-sight velocity distributions (LOSVDs), as it identifies 
the kinematic tracer population ( §5.41) . 



2.1.1. Reduction of the J, H, and K-Band Data 

The 2MASS J, H, and if-band images (Figure [IJ top row) were kindly provided to us 
by Tom Jarrett ( Jarrett et al. 2003 ). and extend out to R ~ 8(6 from the centeiQ. For the 



morphological study of the large-scale properties of the galaxy (§3J), the signal-to-noise ratio 
(S/N) was boosted by constructing a composite J+H+K image, which was then smoothed 
with a Gaussian point-spread function (PSF) having full width at half maximum (FWHM) 
of 7" for r < 60" and FWHM = 11" for r > 60". The 2MASS K-band image was reduced in 
a manner similar to that of the J+H+K image. It proved to have sufficiently high S/N to 
give a reasonably adequate measurement of the outer disk. A comparison of the K-band and 
J+H+K-b&nd photometry provides an estimate of the uncertainties; these are dominated 
by the effects of patchy star formation and other irregularities, not by photon statistics. 



1 Throughout the paper, R refers to distance from the center of the galaxy, while r refers to distance along 
the isophotal semimajor axis. 
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Therefore, we decided to discard the J+H+K-b&nd profile as a tracer of stellar mass in the 
dynamical modeling, and to use the K-band profile instead, as being less affected by these 
irregularities. 

The KPNO 2.f-m telescope i^-band image extends out to R « 3' from the center 
(Figure [H bottom row). Comparison of the 2MASS K-band profile with the KPNO K-band 
profile, which has a higher spatial resolution, showed that the 2MASS profile is affected by 
smoothing only at r < 32", so the 2MASS profile was truncated inside this radius. Thus, 
the KPNO profile was used for 2" < r < 32" and the 2MASS profile at r > 32". 



Chary et all fl2000h obtained HST NICMOS images of NGC 4258 in bandpasses that 
approximate J, H, and K. Eric Becklin and Ranga Chary kindly made the reduced im- 
ages available to us (Figur e El), along w i th th e corresponding PSF images. Details of the 
reductions can be found in Icharv et aD feoooh . The HST WFPC2 V-btrnd image f 32X21 
shows extensive dust at r > 2" on the SW side of the center and a few thinner dust patches 
on the NE side. There is also a scatter of faint sources in the bulge, apparently associated 
with some ongoing star formation. In the NICMOS i^-band image, these faint sources are 
not visible, and the dust patches on the NE side are essentially invisible as well. However, 
the dust on the SW side, although less conspicuous in K than in V, remains a substantial 
problem that we had to correct. After considering a number of alternatives, we chose the 
strategy of creating a symmetrized image of the galaxy, whereby the SW side was replaced 
with the NE side flipped (not rotated) across the major axis. As a check, the profiles from 
both the symmetrized and the unsymmetrized K-b&nd image were compared and found to 
agree extremely well at r < 2". Therefore, this symmetrization is not a factor in the V — K 
color gradient discussed in §4.11 We applied the same procedure to the J- and if-band 
images, as well. 



The nucleus of NGC 4258 emits broad permitted lines (IHo et al.lll997), and is classified 



as a Type 1.9 Seyfert. The AGN at the center is very red (I Chary et al 



20001 ) and, in the 



K band, it appears as a bright, unresolved point source that clobbers the central brightness 
distribution of the stars and needs to be subtracted. Unfortunately, the PSF provided by 
Chary et al.l (120001 ) was determined using the TinyTim routine and not from actual images 
taken at the same time as the galaxy observations. It proves to be a mediocre match to 
the AGN at r < 0"2, in the sense that there is a residual diffraction pattern no matter 
how the PSF is scaled. We applied the follow ing procedure to improve the qu ality of the 
PSF subtraction (cf. Figure [3]). Using VISTA (ILauer. Stover, fc Terndruplll983l ). we scaled 
the TinyTim PSF image by various factors and measured the resulting surface brightness 
profiles in the residual image using the unsymmetrized galaxy image. We also examined the 
residual images to see which ones had the small-scale structure, illustrated in the top left 
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image of Figure [31 optimally subtracted. Most of the weight in the choice of the final PSF 
scale factor came from how well the first diffraction ring of the AGN was subtracted. The 
scale factor that we finally adopted in this way yields ellipticity and position angle profiles 
that are both smooth and well behaved, which is reassuring. This is not the case with scale 
factors that differ by more than ±10% from the adopted value. However, ultimately no value 
of the scale factor allows both a convincing subtraction of the first diffraction ring and a 
smooth profile at r < 0"2 that follows the power law seen at larger radii. This is probably 
a result of the imperfect match of the TinyTim-generated PSF to the K-b&nd data. It is 
unlikely that a real galaxy feature is involved, considering that no such feature is detectable 
in the higher- resolution \^-band image in the radial range 0"05 < r < 0"2 (cf. §2.1.21) . 



We also applied forty iterations of Lucy deconvolution (jLucyl 11974: iRichardson 11972) 
to the HST NICMOS J, H, and J\~-band images using the PSFs from Icharv et al.l feoool ). 
Because the PSFs are not ideal, we applied the deconvolution to images with the PSF of 
the AGN subtracted as explained above. Still, the deconvolved K-b&nd image showed bad 
"ringing," so we truncated it inside r = 0"46. The J-band profile appeared better behaved, 
at least partly since the AGN is much dimmer at this bluer wavelength, although there was 
also some sign of a residual AGN at r < O'.'l. The J, H, and K-band deconvolved profiles 
agree very well at r > 0"5, where the AGN is not a problem, and the ellipticity profiles are 
all well behaved. 

Based on the preceding analysis, we decided (i) to use the original (not deconvolved) 
if -band photometry for r > 2"8 (KPNO data at 2"8 < r < 32", and 2MASS data at 
32" < r < 560"), (ii) to average the deconvolved NICMOS K-band and J-band photometry 
for 0V2 < r < 2'/8, and (iii) to use deconvolved J-band data for 0"038 < r < 0"2 in order to 
avoid the AGN as much as possible. The resulting profile is shown in Figure HI However, 
even in J, the points at r < 0"11 are probably still affected by the AGN. We elaborate on 
this point in §4.21 



2.1.2. Reduction of the V and I -Band Data 

We used the MDM 1.3-m telescope to obtain an image of the galaxy in the R band 
(Figure [5]). The image was reduced and dust-clipped in the usual way. The luminosity 
profile was then extracted and matched to the V^-band HST profile, discussed below, in the 
range 5" < r < 7", using a \^-band zeropoint. There were no significant color, ellipticity, 
or position angle (PA) gradients between the \^-band and i?-band profiles over the radial 
range in which they were matched. Henceforth, we will refer to this composite profile, made 
by stitching together the HST V-band profile with the i?-band profile shifted by a constant 
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offset, as the \^-band profile. It extends out to r ~ 150". 



We obtained HST WFPC2 images of the NGC 4258 nucleus through the F547M filter 
under program GO-8591. The F547M filter was selected to exclude any bright emission 
lines associated with the AGN; it is roughly equivalent to a narrow V^-band filter, and is 
referred to as such in the remainder of this paper. The galaxy was centered in the PCI chip 
and dithered in a 2 x 2 square raster of ~0.5 pixel steps over four separate 400 s exposures 
to maximize the spatial resolution of the complete data set. The four dither ed images wer e 
combined into a single Nyquist-sampled "super image" using the algorithm of lLauerl (jl999al ). 
which provides for the optimal combination of undersampled images without any associated 
degradation of the resolution. This final image has a pixel scale of 0'/0228, twice as fine as 
that of the original images. However, the Nyquist image remains modulated by th e HST 
PSF ; thus, it wa s deco nvolved prior to further analysis using 80 iterations of the iLucyi (119741 ) 
and iRichardsonl (jl972j) algorithm. The PSF was estimated using the TinyTim package, but 
incorporates the WFPC2 pixel-response function recovered by iLauerl (jl999bl ). The final 
deconvolved K-band image of NGC 4258 is shown in Figure [6j 

The K-band surface brightness profile of NGC 4258, shown in Figure HI was estimated 
fr om the F547M Nyquist image using a combination of the high sp atial resolution algorithm 
of ILauerl (119851 ) and the least-squares estimator of ILauerl (119861 ). The patchy dust seen 
throughout the bulge of NGC 4258 complicates measurement of the brightness profile; as 
is visible in Figure El there is a compact dust feature just slightly offset from the nuclear 
source that was especially problematic. For all dust features other than the central patch, 
pixels affected by dust could simply be excluded from the least-squares profile estimator. 
In practice, this was done iteratively by comparing the image to a model reconstructed 
from the profile to isolate faint dust absorption that may have been missed during the initial 
inspection of the image. Unfortunately, pixels could not be masked out of the high-resolution 
profile algorithm, which is used in preference to the least-squares estimator for r < 0'/5, so 
the nuclear dust patch was first filled in with an estimated correction provided by the lower- 
resolution least-squares algorithm. This initial correction was then refined by replacing it 
with an improved estimate provided by a model reconstructed by the high-resolution profile. 

Figure M shows the central portion of the l/-band image divided by a model reconstructed 
from the final brightness profile. Apart from the nuclear dust patch, the residuals are flat, 
showing that the brightness profile extracted is faithful to the central structure of NGC 
4258. The nuclear dust patch, itself, is clearly compact, affecting only a single quadrant of 
the nucleus, which again allowed it to be isolated from the profile estimation. As a "sanity 
check" on the brightness profile, we compared the profile to an intensity trace taken along 
the bulge semimajor axis on the side of the nucleus opposite the nuclear dust patch. The 
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two measures agreed extremely well, showing that the final profile is not strongly affected 
by any residual dust features or the pixels excluded from the analysis. 

In addition to the H ST F547M data we also analyzed two archival PCI F791W (J 
band) images obtained by lCecil et al.l (120001 ) under program GO-6563. Unfortunately, both 
images were saturated at r < 0'/09 from the nucleus, and thus provide information only at 
somewhat larger radii. The images were combined and deconvolved; however, without a full 
accounting of the light contributed by the nucleus, the /-band profile should be regarded 
with caution for r < 0'.'5. 

Using the F547M and the F791W images, we created a F547M/F791W color map of 
the nuclear region of the galaxy, approximating V — I (Figure [7]). The color map shows 
more clearly the distribution of dust in the center, and will be further discussed in §4.11 A 
larger-scale view of the central region of the galaxy is shown in FigureEl which is a composite 
color image created by combining the F54 7M and the F791W images along with an archival 
F300W ([/-band) image, also obtained bv lCecil et all j200oh . 



2.2. Spectroscopy 

We used the Ca II triplet absorption line (8498 A, 8542 A, 8662 A; hereafter CaT) to 
determine the stellar LOSVDs at the center of NGC 4258, and at various positions along 
the major axis (out to R = 18'.'18) as well as along the minor axis (out to R = ll'/69). The 
spectra were obtained from the ground as well as with the HST Space Telescope Imaging 
Spectrograph (STIS). Spectrograph configurations are shown in Table [2J 

In particular, we made 24 cosmic-ray-split exposures of NGC 4258 with STIS to obtain 
high-quality nucleus-centered CaT spectra along the major axis of the galaxy, out to R = 
0'/80. The STIS slit was placed at PA = 140° instead of directly on the major axis (PA = 
150°) because of gui de-star constraints. Our reduction procedure closely followed that of 



Pinkney et al.l (120031 ). The procedure bypasses the pipeline ("Calstis") reduction so that we 
can employ "self-darks," i.e., dark frames that are constructed by combining our dithered 
datasets with the galaxy light rejected. Our method works best if all of the exposures are 
taken within a short period of time because of the rapidly changing dark pattern on the 
STIS CCD. This was not quite the case here, since the exposures were taken on March 12 
and 21, 2001. Nevertheless, a single self-dark worked well. 

There is one s mall d ifference between our STIS reduction procedure and that described 



by iPinkney et al.l (120031 ): the initial dark subtraction was done with an archival, "weekly" 



dark rather than by a sum of the galaxy frames with a region around the galaxy masked out. 
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This created a smoother dark on which to begin further iterations using the same "self-dark" 
strategy. 

We also used the Modspec spectrograph on the MDM 2.4-m telescope to obtain ground- 
based CaT spectra along both the major and minor axes at larger radii. The configuration 
is again shown in Table [21 The seeing varied between 1" and l'/5 (FWHM). Kinematic 
parameters were extracted out to R = 18'.'18 and R = 11'.'69 along the major and minor axis, 



respec tively. Again, the reduction procedure followed closely that described by lPinkney et al. 



(|2003|) 



Figures [9] and [10] show spectra extracted at several radial positions from the STIS and 
the ground-based data, respectively. The emission-line feature visible near 8620 A (rest 
frame) in the nuclear STIS spectrum in Figure is most probably the Fe II 8618 A line, 
and had to be removed before the LOSVD fitting could be done. Consequently, there was a 
concern that the central (or central few) LOSVDs would be less reliable, but it proved easy 
to interpolate under the line (see Figured]), so in the end this was not a problem. 

We decided to discard the central ground-based LOSVD along the major axis, because 
of concerns about template mismatch and AGN contamination. In retrospect, the central 
ground-based LOSVD on both the major and minor axes should have been dropped since 
the STIS data had comparable statistical quality and much better understood PSFs. In fact, 
we used the minor axis central LOSVD, except in one experiment, where its exclusion had 
no effect (see §5.3p . 

In a dust-free axisymmetric galaxy that also exhibits symmetry about the equatorial 
plane, such as assumed by our modeling method, kinematic quantities also manifest symme- 
tries about the center. This entitles us to "symmetrize," with respect to the center of the 
galaxy, the LOSVDs along the major and minor axes, for both STIS and Mo dspec. This pro 



cess h elps alleviate potential biases in the extraction of the velocity profile. iGebhardt et al. 



(120031 ) provide more details on how this was done. Unfortunately, NGC 4258 contains sub- 
stantial quantities of dust, especially near the dynamically important central region (Figures 
El E]), which calls into question the validity of this assumption. Potential implications are 
discussed in 361 



We used the maximum penalized likelihood method (MPL), as described by lPinkney et al 



(120031 ). to derive LOSVDs from all the CaT spectra. A nonparametric form of the LOSVDs, 
rather than Gauss-Hermite moments, was used for the dynamical modeling. Each LOSVD 
was represented by 13 equally spaced velocity bins, with the uncertainty in each bin deter- 
mined from Monte Carlo simulations. 



Figure fTTl illustrates V, a, if 3, and if 4, the fitting parameters from the truncated Gauss- 
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der Marel fc Franxll993L 



Hermite series expansions of the LOSVDs (IGerhard et al.lll993l ; lvan 
along the major and minor axes of the galaxy, both from STIS and from the ground. The 
same parameters in tabular form can be found in Table [3J We stress again that a nonpara- 
metric form of the LOSVDs was actually used for the dynamical modeling. We anticipate 
our best fitting model discussed in §5.1 by illustrating its projected Gauss-Hermite moments 
here. 



3. THE LARGE-SCALE MORPHOLOGY 



NGC 4258 is a typical oval-disk galaxy (jKormendylll982l ). classified as SAB(s)bc. Mor 



phologically, the disk consists of two nested regions of slowly varying surface brightness: the 
"inner oval" component (PA ~ 156°) is much brighter than the "outer oval" component (PA 
~ 150°), and both have relatively sharp outer edges, at R ~ 200" and R ~ 560", respectively 
(Figured]). Embedded in the inner oval component is a bulge (PA ~ 146°) which remains 
important out to r ~ 40", as well as a weak bar (PA ~ 11°) extending out to R ~ 150". 
These components can also be identified using the profiles in Figure HI 

We are interested in the large-scale morphological properties of NGC 4258 primarily to 
establish that there are no major deviations from the assumption of axial symmetry, which 
is built into our modeling method ( §5.11) . Deviations from axisymmetry can be important for 
two reasons. First, orbits in a non-axisymmetric potential do not conserve any component 
of the angular momentum. Therefore, the orbital structure of a non-axisymmetric system 
can be considerably different from that of an axisymmetric one. Second, non-axisymmetric 
deviations can induce kinematic signatures in the radial velocity profile similar to those 
created by a central mass concentration. Either effect can bias a black hole mass estimate 
made under the assumption of axisymmetry. 

The ellipticity profiles in V and in the NIR are quite different in the innermost 0"2, but 
they agree rather well further out. We consider the ellipticity profile in V to be overall more 
reliable than that in the NIR due to the higher resolution in V. It becomes progressively 
flatter from the center out to r ~ 5", and then it remains approximately constant between 
r ~ 5" and 12". The ellipticity starts dropping again out to r ~ 40", but this is probably 
caused by the weak bar, which adds extra light in the minor- axis direction while probably 
having no effect on the major-axis profile. The onset of the dominance of the disk is signified 
by the increasing ellipticity outward of r ~ 40". 



If we wanted to rigorously test these morphological characterizations, we would have to 
perform a full bulge-bar-disk photometric decomposition. However, such a decomposition 
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would be very uncertain because of parameter coupling: the bar is very weak, and trading 
light between the components while correspondingly changing their isophote shapes would 
provide considerable freedom to tinker with the bulge ellipticity. 

In the absence of extinction, the projected isophotes of a galaxy that has emissivity 
constant on coaxial and similar spher oids (as assumed in our dynam ical modeling) are con- 
centric, coaxial, and similar ellipses (jContopoulosI Il956l ; iFishl Il96ll ). The absence of any 
strong isophote twists out to r ~ 40", the beginning of the weak bar, indicates that, to a 
good approximation, NGC 4258 is axisymmetric out to that radius. Furthermore, the fact 
that the PA of the bulge and of the outer disk agree quite well is evidence against a triaxial 
bulge. However, there does exist a strong isophote twist between 40" and 150", caused by the 
bar, while Figure [TT] shows some minor-axis rotation in the bulge, outside R ~ 10", which 
could be indicative of triaxiality at intermediate radii. 

Even though our current algorithm cannot model triaxial components, the dynamical 
effect of the bar on the measured value of the black hole mass is probably small, for at 
least three reasons: (i) M. is mostly determined from the dynamics near the center, where 
triaxial deviations are least significant; (ii) bars generally have relatively little effect on the 
kinematics — for example, even t he much stronger bar o f NGC 936 has little effect on the 
velocity dispersion of that galaxy (jKormendylll983l . Il984l ); (iii) and finally the velocities of 
bar stars that cross the slit, while probably somewhat bigger than the local circular velocity 
(a ~5-10% effect), are seen moving nearly perpendicular to the line-of-sight because the bar 
is strongly inclined to the major axis, and hence they probably do not signicantly affect the 
observed kinematics. 

In summary, NGC 4258 contains non-axisymmetric structures which probably do affect 
the kinematics, and hence the black hole mass estimate, but probably only in a minor way. 



4. VARIATIONS IN M/L RATIO 

Stellar population gradients, dust extinction, and contamination by an AGN can all 
induce variations in the M/L ratio, which have to be taken into account when deriving the 
stellar mass density profile from the luminosity density profile. We first discuss variations in 
M/L far from the AGN-dominated nuclear region, and then we describe the complications 
due to the AGN. 
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4.1. The Stellar M/L Ratio Far from the AGN Point Source 

Color gradients in the outer disk of NGC 4258 are small, even between the R and K 
bands (Figure Hj), implying that the M/L ratio in the disk is not seriously affected by star 
formation or dust. However, inside ~2" from the nucleus, the l/-band profile "peels off" from 
the NIR profile and becomes flatter. This color gradient can be seen in Figure HI but becomes 
more apparent in Figure [T2l which shows the inner 10" of the surface brightness profile. Also, 
the color map shown in Figure [7] clearly identifies the color gradient as a "red," diffuse area 
surrounding the nucleus out to r « 2" (R m 1" — 2"). This is an issue for the dynamical 
modeling because it takes place in the dynamically important central region, which affects 
the entire model: we must determine which, if either, of the two profiles (V-band or NIR) 
traces stellar mass. 

Could the color gradient be caused by variations in the stellar population? The strongest 
evidence against variations in the stellar population is the absence of any marked color 
gradient between the J and K profiles, which should be present if the central NIR brightening 
were caused by either a metallicity gradient or red supergiant stars. Furthermore, although 
some traces of scattered star formation can be seen on the l^-band image (Figure Ej) as well 
as on the color map (Figure [Tj), they are absent in the NIR images (Figure [2]), even after 
accounting for the different FWHM of the PSFs. 

More information can be obtained from the equivalent width (EW) profiles of the three 
CaT line components, which we computed from the STIS spectra extending out to l'/5 
from the center (Figure [TBI . A decreasing EW profile would be evidence of a young stellar 
population or of a brightening /-band nonthermal continuum. A drop is seen in the central 
(X'3, and is discussed in the next subsection. Here we are concerned with possible trends 
beyond that radius. For r > 0'.'3, there is a hint of a slow increase in the EW away from 
the center, but the data are also consistent with a flat profile. However, the datapoints at 
r > 0'.'8 suffer from the low S/N of the originating spectra, and we have no data as far out 
as ~2", where the F-band profile begins to peels off. 

The other possible explanation is the presence of diffuse dust, which would produce a 
small gradient m J — K but a larger one in V — J, just as observed. The observed smoothness 
in the V — J color gradient would then imply a similarly smooth increase in the projected 
density of diffuse dust toward the center. 

In the end, we were unable to determine unequivocally whether the origin of the color 
gradient is a variation in the stellar population or the presence of a diffuse dust compo- 
nent. We decided to construct models for both of these possibilities, and investigate the 
uncertainties induced on the black hole mass ( §5.4p . 
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4.2. The Stellar M/L Ratio Near the AGN Point Source 

Dynamical modeling only strongly constrains the total mass, M to t(ro) = M. + M*(ro), 
inside a small radius ro around the center, where M*(ro) is the enclosed stellar mass, and ro 
is comparable to the spatial resolution of the nuclear spectrum (~0 / .'l). Therefore, there is a 
degeneracy between M. and M^(r ), which can be resolved only if we know M*(ro) separately, 
i.e., only if we know the luminosity profile and the stellar M/L ratio inside ro- The latter 
is assumed to be constant throughout the galaxy and is predominantly constrained by the 
dynamical model at large radii ( §5.11) . In essence, M. = M tot (r ) — M*(r ) corresponds to 
the "excess" dynamical mass inside r after subtraction of the mass due to stars. Therefore, 
we must again determine which, if either, of the two luminosity profiles (V band or NIR) 
traces stellar mass, M*. 

The l/-band profile, which is well determined down to r ~ O'/Oll, manifests a progres- 
sively more shallow slope for r < 0'.'15, interrupted by a peak within 0'.'04 of the nucleus. 
Although the cause of the peak is not entirely clear, it is unlikely that it is due to light from 
an old stellar population. The CaT EW profile (Figure [13]) shows a ~30% decline within 
Of! 2 of the nucleus, which strongly suggests contamination by a continuum source, possibly 
OB stars or a nonthermal AGN. 

Notwithstanding the uncertain origin of the peak, it seems unlikely that it traces stellar 
mass. Therefore, we replaced the innermost three l/-band datapoints in a way that extrap- 
olates inward the profile at 0'.'04 < r < 0'/15 (Figure [121 §5.4p . However, we also created 
models with the central peak, to investigate its effect on M.. Not surprisingly, the effect was 
negligible because the light in the peak corresponds to a very small mass ( §5.4p . 

In §2.1.11 we saw that, at r < O'.'ll, the J-band profile begins to show signs of contam- 
ination by AGN light, and the PSF of the AGN could not be subtracted reliably. We also 
saw earlier that the central depression in the CaT EW profile begins at r ~ 0'.'2, and implies 
a brightening of the CaT continuum, which lies in the I band, inside that radius. The AGN 
is red, so it seems likely that the J-band profile is even more affected by the AGN continuum 
radiation than the /-band profile. While luminous AGNs normally have a blue featureless 
continuum, low-luminosit y AGNs such as NGC 4258 tend to have emission spectra that 



peak in t he mid-infrared (lHolll999l ). probably as a consequence of their low accretion rates 
(IHol 120081 ) . Taking all these considerations into account, we decided that the J-band profile 
inward of r ~ 0'/2 cannot be trusted to trace mass (i.e., the old stellar population). Accord- 
ingly, inside r ~ Of! 2 we replaced it with a profile that follows the V^-band profile without the 
central peak. As explained above, this profile is presumably not affected by the AGN. The 
related dynamical models are discussed in §5.4[ 
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5. DYNAMICAL MODELING 
5.1. The Method 

A detailed account of our dynamical modeling method is provided in the Appendix. 
Here we present a short description with emphasis on NGC 4258. 

The first step involves a deprojection of the surface photometry to obtain the three- 
dimensional luminosity density. This is done under the assumption that emissivity is constant 
on similar, c oaxial spheroids . Th e deprojection is performed via a nonparametric Abel 



inversion (cf. iGebhardt et al.l Il996l ) assuming a value for the inclination angle (i = 72° in 
the base model). The gravitational potential can then be recovered by specifying a central 
point mass ("black hole") M. and the stellar M/L ratio, T, assumed constant throughout 
the galaxy. The dynamical models for NGC 4258 extend out to r = 100". 

The dynamical models are constructed for a grid of (M., T) values. Each model consists 
of a superposition of representative orbits, appropriately weighted to match the observed 
kinematics, subject to constraints that force the stellar density distribution to match the 
luminous density distribution of the galaxy. Each model is computed using iV w 7000 
orbits. Orbits are integrated for about 100 crossing times. 

Self-consistency is enforced by requiring that the model satisfy exactly the photometric 
constraints on a polar grid of 15 radial shells as described in the appendix, each subdivided 
into 4 bins in polar angle. The grid used to solve Poisson's equation is 4X finer than the 
grid used to bin the kinematic results. For most orbits energy is conserved to much better 
than 1%. 

For each model, the orbit weights that best fit the observed kinematics are determined 
by minimizing 

X 2 = 5^ [ lk >° ~ S m ik W i l a l > W 
k 

where lk )0 is the light at a specific position and velocity (the index k runs over both position 
and velocity) with uncertainty a^, and where is the mass deposited by the i th orbit 
weighted by Wi in the k th projected position and velocity bin; cf. §2.21) . We use a maximum- 
entropy technique as described in the Appendix to minimize \ 2 using non-negative orbit 
occupation numbers. 

For each dataset, we determine x 2 (M,,T). The minimum determines the best (most 
probable) combination of M. and T for the galaxy. The A% 2 = 1 contour on the (M.,T) 
plane, where A% 2 = \ 2 ~ Xmim determines the nominal "la" uncertainty ranges for M. and 
T. 
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The process of converting observational quantities to model input entails assumptions 
and systematic uncertainties that have to be taken into account when determining the un- 
certainty of the "best" M., which may be larger than the nominal "la" confidence band on 
the (M., T) plane. Therefore, the previous steps must be repeated for a number of plausible 
parameter choices (other than M. and T) in order to determine the overall uncertainty. We 
cannot afford to make an exhaustive exploration of the parameter space because of the time 
required to create the models. Instead, we opt to construct a base model corresponding to 
our "best guess" for the parameter values and investigate the effect on M. of varying each 
parameter independently. 

The results from the dynamical modeling are summarized in Table HI and the corre- 
sponding x 2 contour maps are shown in Figure [T41 The black hole mass for the base model, 
which uses the l/-band profile without the central peak for both the mass and the kinematic 
tracer profiles, is M. = 3.3lj^i7X 1O 7 M with a total (not per degree of freedom) Xmm ~ 290. 

The total number of parameters is 27 (LOSVDs) t imes 13 (velocity bins) , which is 351. 
Typically, the reduced x 2 is 0.3 to 0.5 in our models (IGebhardt et al.l 120031 ) . which is less 
than the expected value of unity because of covariance in the LOSVD components. For NGC 
4258, the reduced x 2 is near unity, but most of the contribution to x 2 comes from the minor 
axis. The velocity centroid on the minor axis varies much more than the stated uncertainties 
(see Figure 11), whereas an axisymmetric model would have zero velocity along the minor 
axis. Regardless of whether this contribution is due to underestimated uncertainties or real 
non-axisymmetric motion, it is the main cause for the inflated x 2 ■ Along the major axis, the 
reduced x 2 is about 0.4, which is typical for the orbit-based models. 



5.2. Ambiguities in the Deprojection Parameters 

5.2.1. Isophotal Ellipticity (e) 

The (projected) ellipticity profile of NGC 4258 is discussed in §|3j The profile varies with 
radius, but the origin of the variation (disk, weak bar, genuine variations in the axis ratios of 
the bulge isodensity spheroids) is not clear. Therefore, although it is possible to incorporate 
an arbitrary ellipticity profile e(r) in our modeling algorithm, we preferred instead to make 
models with a constant ellipticity profile, and examine the extent to which M. is affected by 
the adopted value. 

For the base model we used e = 0.35, corresponding to the mean ellipticity in the radial 
range 2" < r < 65". The effect of varying e was examined in model Dl, for which we used 
e = 0.45. We obtain M. = 3.48lg;3g x 1O 7 M , a slightly higher mass but within the same 



"la" confidence as the base model. 



Nonetheless, the assumption of constant ellipticity is clearly problematic because the 
ellipticity profile varies between ~0.1 and ~0.6. 



5.2.2. Disk Inclination Angle (i) 

We need to know the inclination of the disk to deduce the intrinsic ellipticity of the 
isophotes from their projected ellipticity, under the assumption that the equatorial plane of 
the bulge of the galaxy is coplanar with the galactic disk. 

The outer disk (beyond r m 200") has average and maximum ellipticities of e = 0.62 
and 0.64, respectively. If this disk is infinitesimally thin, the corresponding inclination 
angles are i = 69° and 70°, respectively, where edge-on corresponds to i = 90° . If this 



disk were thick with an intrinsic axis ratio of 0.25 (ISandage. Freeman, fc Stokeslll97Cll ). then 



the corresponding i nclinations would b e i = 73° and 75°, respectively. At large radii, the 



optical PA is 150°. [Van Albadal (119801 ). using H I kinematic observations, derives the same 
PA and i = 72°, which is consistent with the photometrically derived inclination given the 
uncertainty in the intrinsic thickness of the disk. 

We therefore adopted i = 72° for the base model, but we also examined the sensitivity 
of M. to variations in i. Model D2 is identical to the base model except i = 62°, a value 
close to the lowest inclination angles that we found in the literature. We obtained a ~10% 
higher mass, M. = 3.62 jlo'll x 10 7 M Q which is consistent with a M. oc (sinz) -1 dependence. 
This would be expected in a rotationally supported system. It appears from our models that 
the inner part of the galaxy is rotating quite rapidly. For this model, Xmm = 303, which is 
greater than 280, the Xmin f° r the base model. This could be interpreted as evidence that 
the base model is more likely to be true, at least within the framework of assumptions of 
the modeling method. 

We did not attempt to obtain a better estimate for the errors due to deprojection (by 
running models for more combinations or % and e) because the uncertainties in the luminosity 
profiles ( §5.4p dominate the errors. 



5.3. Ambiguities in the Kinematic Data 



In §2.21 we mentioned that the presence of an emission line near the CaT in the nuclear 
STIS spectrum may have affected the reliability of the LOSVD extraction. In fact, the next 
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two STIS spectra (at R = (X'05 and O'.'IO) show (weaker) signatures of the same emission 
line, as well. We also mentioned that we have some concerns about the quality of the central 
LOSVD along the minor axis, in the ground-based data. In this section we remove the 
suspect LOSVDs from the x 2 fit, and examine the effect on M.. 

In models Kl and K2 we isolate the uncertainty introduced to M. by the potentially 
unreliable LOSVD extraction. Model Kl is identical to the base model, except the nuclear 
STIS LOSVD is not included as a constraint. In model K2, the three innermost STIS 
LOSVDs are not included. In model Kl, the effect on M. is very small. This is reassuring, 
since it is the nuclear spectrum that would have been affected the most. Model K2 yields 
a significantly lower M. = 2.20lg'3i x 10 7 M q but with a much higher uncertainty, which 
overlaps with the base model at the ~1.5<r level. 

An interesting question is whether HST spectroscopy is required to find evidence for 
a black hole. Using ground-based kinematic data alone should, of course, yield a black 
hole mass consistent with the mass obtained by also including HST data, albeit with larger 
uncertainties, depending primarily on the angular size of the radius of influence of the black 
hole on the sky. We address this question for NGC 4258 with model K3, which is identical 
to the base model, except only ground-based kinematic data are used. It turns out that the 
stellar M/L ratio can be constrained quite well, but the black hole mass, M. = 1.03^Q28 x 
1O 7 M , is very uncertain. This is not surprising, considering that the radius of influence for 
the black hole in the base model is GM./<7g ~ (X'32, where a e = 105 km s -1 is the velocity 
dispersion in the main body of the bulge. This is substantially smaller than the spatial 
resolution of the ground-based kinematic data, for which FWHM = 1"— 1'/5. Therefore, HST 
kinematic data are critical for a well-constrained estimate of the black hole mass in this 
galaxy. HST photometric data are also critical for galaxies where the luminosity profile near 
the center cannot be reliably extrapolated from the profile farther out, as is also the case 
with this galaxy (unfortunately, it is still not possible to determine unambiguously the inner 
mass profile, as discussed in §5.41) . 

We also constructed model K4, which uses all kinematic data (STIS+Modspec) except 
the central ground-based LOSVD along the minor axis (the central ground-based LOSVD 
along the major axis was not included in any model). We find that M. is essentially identical 
to that of the base model. 
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5.4. Ambiguities in the Mass and Kinematic Tracer Profiles 

In §4] we saw that stellar population gradients and/or dust extinction in the innermost 
~2" of the galaxy, as well as light from the unresolved AGN, can induce variations in the 
measured stellar Mj L ratio near the center of the galaxy. These have to be taken into account 
to determine which, if any, of the available luminosity profiles traces the mass density profile, 
and which, if any, corresponds to the luminosity profile of the kinematic tracer population. 
The latter profile is needed for the luminosity weighting (normalization) of the LOSVDs. 
Unfortunately, we were unable to make these determinations convincingly, for the reasons 
explained in §HJ 

Consequently, we decided to identify three possible luminosity profiles, and to make 
models that use them as either mass or kinematic tracer profiles. The expectation is that 
these three profiles correspond to limiting cases, and that the true profiles lie somewhere in 
between. Then, the spread in M. that we obtain from these limiting cases would be a good 
indication of the M. uncertainty caused by our inability to determine the true profiles. 

The profiles, which are depicted in Figure [12j are the following: 

(1) The V profile: This is identical to the U-band profile as measured. 

(2) The V corr profile: This is identical to the V profile, with the three innermost points 
replaced by a shallower core, as explained in §4.21 

(3) The NIR COTT profile: For r > (X'2, this is identical to the spliced J, J+K, and 
K profiles in the top panel of Figure HI Inward of (X'2, that profile is replaced by one that 
parallels (shows no color gradient with respect to) the V corr profile, for the reasons mentioned 
in §E2 

If the observed color gradient is mostly or entirely due to the presence of diffuse dust or 
metallicity, then the NIR corr profile is a better tracer of mass than the V or the V COTT profiles. 
However, the l/-band profile has higher intrinsic resolution and a dimmer AGN, and the 
WFPC2 is better understood photometrically than NICMOS. For these reasons, we try all 
three versions as stellar mass profiles. 

The kinematic tracer population was observed in the CaT wavelength, which lies in the 
/ band. It would thus be natural to use the /-band profile as the kinematic tracer profile. 
However, we saw in §2. 1.21 that the J-band image is saturated in the central 0'/09, and thus 
cannot be deconvolved. The NIR profile shares the same slope with the /-band profile at 
larger radii (where / is not saturated), and could arguably be used instead. Unfortunately, 
it is also ill-determined near the center due to the AGN and reduced resolution ( §4.ip . 
Consequently, we decided to again test both NIR corr and V COTT as the kinematic tracer profile, 
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and examine the induced uncertainty in M.. For the base model we adopted V COII because of 
its higher resolution and smaller AGN contamination. This is somewhat arbitrary, but the 
choice of the kinematic tracer profile turns out to have a very small effect on M., as we will 
see below. 

These scenarios were tested in models PI through P4. Models P1-P3 substitute NIR corr 
for V COTT in the mass and/or kinematic tracer profiles, and show that the choice of the stellar 
mass profile matters much more than the choice of the tracer profile. When V COTr is used for 
the mass profile (base model and P2), M. = 3.31+°;^ x 1O 7 M or M. = 3.33^ x 1O 7 M , 
depending on whether V COTT or NIR corr is used, respectively, for the tracer profile. These 
numbers become 2.49l£™ x 1O 7 M and 2.43l{jJJ x 1O 7 M , respectively, when NIR corr is 
used for the mass profile (models PI and P3). The relative insensitivity to the choice of 
the tracer profile is not surprising, considering that, e.g., for a spherical system, the tracer 
profile v{r) affects mass M(r) within radius r only via dlnu/dlnr, which by construction is 
nearly identical for both V corr and NIR corr near the center of the galaxy. 

Model P4 substitutes V for V covv in the mass profile. It yields M. = 3.25t° H x 1O 7 M , 
very similar to the base model, reflecting the fact that the only difference between P4 and the 
base model is the small luminosity excess within 0'.'04 of the nucleus, well below the spatial 
resolution of the kinematic data (~0 / .'l). The stellar M/L ratio is constrained at larger radii, 
so that the only difference in the stellar mass profile, M*(r), between P4 and the base model 
is due to this small luminosity excess. From the M. — M*(ro) degeneracy, discussed in §4.2^ 
the only option for the modeling algorithm is to reduce M. by a (small) amount equal to the 
difference in stellar mass inside tq px O'/l between P4 and the base model. 

This point is illustrated in Figure [16j which shows the stellar mass profile, M*(r), 
and the total mass profile, M. + M*(r), for the base model and for models PI and P4. 
The mass profiles of these models correspond to the three profile choices (Korr, NIR corr , 
and V, respectively) discussed above. It is clear that dynamic modeling assigns black hole 
masses to P4 and to the base model such that they both contain the same total mass inside 
r ~ 0'/04 (since both models are constrained by the same kinematic data and the same V- 
band mass profile beyond r ~ 0'/04). This is not the case with model PI which, although 
again constrained by the same kinematic dataset, is characterized by a different mass profile 
(NIR corr ) outside of ~ O'/l, the spatial resolution of the kinematic data. 
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6. DISCUSSION 

The determination of the central black hole mass in NGC 4258 from stellar dynamics 
proved harder than anticipated. The main source of difficulty (and uncertainty in M.) comes 
from the presence of a V — J color gradient in the central ~2" of the galaxy, which prevents 
us from determining the stellar mass profile near the center with confidence. The model 
that we trust most ("base model" in Table H]) yields M. = 3.3ll !i7 x 1O 7 M and assumes 
that the stellar mass is traced by the "corrected" l^-band profile, or V CO rr (Figure [T2|) . This 
is identical to the V^-band profile except for the removal of a small luminous "bump" in the 
central ~0'.'04, which we attribute to the AGN. Using the uncorrected y-band profile yields 
a very similar M. = 3.25±°;^ x 10 7 M Q (model P4 in Table H]). This small black hole mass 
"deficit," compared to the base model, corresponds approximately to the mass of stars in the 
luminosity "bump" and is a consequence of the degeneracy between stellar mass and black 
hole mass very near the center ( §4.2( Figure [ToT) . 

Using the steeper J-band (NIR) profile as the mass tracer (model PI), we obtain M. = 
2.49^Q l j ) o x 1O 7 M , a factor of 25% lower than in the base case. Although V-band light is more 
likely to be compromised as an indicator of the stellar mass distribution than J-band light, 
due to extinction and metallicity gradients, we have greater confidence in the base model 
because (i) the V-band profile is better determined near the center owing to the higher 
resolution of the V images, (ii) the AGN is considerably fainter in V than in J, making it 
easier to subtract, and (iii) the minimum x 2 for the base model is 280, significantly less than 
the minimum \ 2 for the PI model using J-band data (which was 303). 

The spread in M. introduced by uncertainties in the deprojection parameters (models 
Dl and D2) is only of order 10%, so it is the uncertainties in the stellar mass profile that 
dominate the errors. 

After all these sources of error are taken into account, our black hole mass determination 
is ~15% lower than the "preferred" maser determination [(3.82 ±0.01) x 10 7 M Q ] and only 7% 
lower than the lowest maser determination [(3.59 ±0.01) x 10 7 M Q ] (rescaled to our distance 
for the galaxy). We view this level of agreement (2a) as evidence that stellar dynamical mass 
determinations using similar methods to those in this paper, and with data of comparable 
quality, are accurate at the 2a level or better. Our previous work has sought to avoid 
galaxies with dusty centers such as NGC 4258. In this case obscuration of high-velocity 
stars by the compact nuclear dust patch ( §2.1.2p . possibly associated with the masing disk, 
could be responsible for an underestimate of the black hole mass in our work. To these 
problems should be added possible systematic effects due to deviations from axisymmetry, 
which cannot be properly treated by our axisymmetric modeling code. We also note that 
the discrepancy between the mass determination from stellar dynamics and the maser mass 
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determination in NGC 4258 is smaller than the discrepancy between the stellar dynamical 
mass determination and the mass deter mination from individual orbits for our own G alaxy, 
which amounts to about a factor of two (jChakrabarty fc Sahall200ll ; iGhez et al.ll2005l ). Most 
concerns about black hole mass estimates from stellar dynamics have stressed the likelihood 
that the masses are overestimated, whereas here the stellar dynamical mass is smaller than 
the maser mass. 



Finally, although we regard the maser mass in NGC 4258 as the gold standard in black 
hole masses, it is conceivable, however unlikely, that some unrecognized effect has led to an 
overestimate of the mass determined from maser kinematics. 
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A. Description of the Orbit-Superposition Method 



For over two decades, orbit-superposition methods have been used to stu dy the dynam- 
ics of galaxies. Originally invented to construct a model of a triaxial galaxy (jSchwarzschild 
19791 ). the methods have become the tool of choice for interpreting data in terms of equi- 
librium galaxy models. The most important feature of orbit-superposition models is that 
they provide a constructive proof that a given set of photometric and kinematic data can 
(or cannot) be reproduced by an equilibrium stellar system in a specified gravitational field. 
The comput er code described in t his appendix has been used to estimate black hole masses 

ative contributions of dark and 
It is a descendant of the 



in galaxies (IGebhardt et al.ll2003l ) and to determine the re 



luminous matter in elliptical ga laxies (Thomas et al 



2004 



20051 ) 



spherical program described in iRichstone &: Tremaind (119881 ). but improved in two major 
respects: it treats the galaxy as axisymmetric rather than spherical, and it matches the full 
LOSVD of the galaxy at specified positions on the sky, rather than the second moment only. 
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Sim ilar programs have b een developed by Ivan der Marel et al.l (119981 ). ICretton et al.l (119991 ) 
and IValluri et al.l (bopj ). 



A number of preliminary steps must be taken to convert the reduced data into inputs 
for the orbit-superposition method. These include the deprojection of the observed surface 
brightness to construct a light distribution in three-dimensional space and the reduction 
of spectra at different locations on the sky to projected velocity distributions. Both of 
these activities require additional assumptions or choices. We normal ly deproject under the 
assumption that level surfaces in stellar density are coaxial spheroids (jGebhardt et al.lll996l ) 
and we normally use a penalized maximum -likelihood estimator to construct the LOSVDs 



( IGebhardt et al.ll2000a 



Pinkney et al.ll2003l ). Next we compute the gravitational field from 



the three-dimensional light distribution under the assumption that the mass consists of a 
black hole of mass M. and stars having a Mj L ratio independent of position. Dark matter 
can be incorporated into such a model, but we have not done so here. 

Once the preliminaries are complete, the calculation of a dynamical model by this 
method consists of two steps. First, a library of orbits is constructed using initial con- 
ditions that cover all possible locations in phase space. Each orbit's contribution to the 
surface brightness and projected velocities is logged. Second, the orbits are combined to 
match the light distribution and LOSVDs of the galaxy as well as possible (using x 2 param- 
eter to assess the goodness of fit). This procedure is repeated for a set of black hole masses, 
and limits on the latter are set from \ 2 as a function of M.. Note that although different 
orbit libraries are required for each different ratio of Mj L to M. , a single library can be used 
for various values of M. so long as M/L scales with M.. 



A.l. The Orbit Library 

We construct the model on a grid in coordinates r, rj where r is the radius and r\ is the 
latitude. We restrict ourselves to potentials (and mass distributions) symmetric about the 
equator. The grid is evenly spaced in v = sin (77) from v — to 1, and in r in even intervals 
in k defined by 

k = -logfl + V). (Al) 
a \ J 

Bins in (r, r]) are bounded by the grid points. Below we also use standard cylindrical coor- 
dinates (to,z). It is essential to explore the phase space of orbits with a resolution at least 
as fine as the resolution used to construct the models. 

Orbits in axisymmetric potentials always have two isolating integrals, the energy E, and 
the z-axis angular momentum J z . Regardless of the possible presence of a third integral, 
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conservation of E and J z confine the orbit within a region of (w, z) space defined by a 
boundary where the orbits' velocities are zero: 

<$>( WiZ ) + Ji-<E. (A2) 
2 w 

The allowed volume intersects the equator only over a limited range in radius with an upper 
and lower limit. We choose (E, J z ) pairs so that there is an orbit with upper and lower limits 
in the middle of each radial bin. We associate an area with each grid point (AE x AJ Z ) 
by halving the distance to adjacent grid points. For each specified (E,J Z ), the section of 
allowed phase space (we ignore because of the conservation of J z ) that lies on the equator is 
a two-dimensional surface with coordinates w and v m (a surface of section). A regular orbit is 
de fined by a curve on th is surface. We systematically tile this surface in the manner invented 



by lThomas et al.l (120041 ) . launching an orbit from the center of each bin and assigning a phase 



space volume to each orbit by 

M-^/n^ltox*. (A3) 



based on the work of lBinney. Gerhard, fc Hutl (119851 ). where T is the time between successive 



crossings of the equatorial plane. 

The procedure above gives a starting position for each orbit in the equatorial plane of 
the model with a well-defined zu, v w , E, and J z . We use E to determine v z , and we integrate 
the motion of the orbit in the (w, z) plane. Assuming that all orbits cross the equator, it is 
a complete survey of orbits. We follow the orbits of stars in spherical coordinates, to take 
advantage of the near sphericity of the mass distribution in the region where a black hole 
may dominate the gravitational field. Using r as the radial coordinate, rj as the equatorial 
angle, and as the axial angle, the equations of motion are 



(A4) 



dv r + V\ d$ 

dt r dr' 

dv v v t)? 1 <9$ 

_p_ _ — — rj_ _ ^ (A.5) 

dt r r r or] 

v$ = (A6) 
r cos 7/ 



The first two equations are integrated using a standard fourth-order Runge-Kutta 
method with variable timestep, and is recovered at each timestep from equation (A6). 
We obtained good energy conservation by letting dt = er/v C i rc (r), where v C i rc is the velocity 
of a circular orbit at that radius, with e between 0.01 and 0.1. Orbits are followed for about 
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100 crossing times. We terminate the calculation of an orbit if its energy error exceeds 1%. 
Typical energy errors are 0.1%. The potential is set equal to zero at oo. 

The accelerations for the orbit integrations are determined by decomposing the assumed 
galaxy density distribution in spherical harmonics. Specializing immediately to axisymmetry 



gives 



with 



^max 



$(x) = -AttG {Air- {1+1) + B ir l ) P(u), 



1=0 



A, 



r r rl 



p(s, v)Pi(v)dv 



p(s, v)Pi{v)dv 



s l s 2 ds, 
s-« +1 h 2 ds, 



(AT) 

(A8) 
(A9) 



and where v = sin//. The gravitational acceleration on the orbit is determined by differen- 
tiating equation IA7I with respect to r and rj. 



We tested our expansion in a Kuzmin-Kutusov model (IDejonghe fc de Zeeuwl Il988l ) 
with a/c = 0.5, which corresponds to a very flat galaxy (about E5). This model has analytic 
density, potential, and accelerations. Truncating the expansion at P 6 accurately reproduces 
the exact accelerations to 1% or better everywhere in the galaxy, so we set l max = 6 in 
equation (A7). 

The additional acceleration from a mass point is straightforward to add to the radial 
component of the accelerations above. We store and model the galaxy in the angular range 
0° < f] < 90°. We store the three lowest-order internal moments in each bin for each orbit 
(mass m, mv, mvv) for later use. 

Simultaneously, the projected line-of-sight velocities are binned in a three-dimensional 
grid in projected radius R, angle (3 (measured from the major axis of the galaxy on the plane 
perpendicular to the line of sight), and velocity v. The spacing in R and (3 is defined as 
above. Since the orbits are being followed in r and rj, <p is not defined. At each timestep, 
in order to project the orbit onto the line of sight, we choose 100 values of 0, evenly spaced 
from to 7r. 



A. 2. Matching the Light Distribution and Observed Dynamics 



The next step is to find the superposition of o rbits consistent with the stel lar mass 
distribution that best matches the observed LOSVDs. Richstone &: Tremainel ( 1988 ) provide 
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an efficient method to match the mass distribution that can be augmented to minimize x 2 - 
Their method maximizes an objective function 

F = -J2 C i( x i) ( A1 °) 

i 

while satisfying a set of constraints 

Yj = y^mjjXj. (All) 



We choose to maximize the objective function 

Norb / \ JVd ji j \2 

i=i x 7 fc=i K 

where Wi is the weight of the i th orbit and Ik is the light in the k th bin in the observational 
phase space composed of the position on the sky and the line of sight velocity. The more 
familiar LOSVD is the distribution of light in line of sight velocity at a single position on 
the sky. Thus the index k identifies a small range in both projected velocity and projected 
position. The variable a is an adjustable parameter discussed further below. There are N or b 
orbits in the orbit library and k varies from 1 to N d , where N d is the number of positions 
at which the LOSVD has been measured times the number of velocity intervals at each 
measurement. The first sum in equation (A12) is an entropy, and the second is a \ 2 that 
measures the goodness of fit of the model LOSVDs to the observed LOSVDs. The sum of 
the orbits must also match the deprojected (3D) light distribution in the N bin bins: 

^orb 

M j = ^2 m i3 W i for 3 = 1 ' N f»ri, (A13) 
i=l 

where m^- is the mass contribution of the i th orbit to the j th bin. 

In order to put equation (A12) in the form of equations (A10) and (All), we make the 
following substitutions and identifications. First, we define Xj = for i — 1, N orb , and add 
an additional set of variables that measure the error in the light projected into each LOSVD 
bin XN orb+ k = Ik ~ h,o for k = l,Nj. Since the light in each bin of the model LOSVDs is 
given by the contribution of each orbit to that projected velocity and position there is an 
additional set of constraint equations of the form the form 

^orb 

h,o = ^ m H W i ~ X 3 for J = N bin + 1, N bin + N d , (A 14) 

i=l 
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where k = j — N orb and rriij is the contribution of the i th orbit to the k th component of the 



set of LOSVDs. 
Then we have 



Q = Xi log (xi/AQi) for % < N orb 

C t = ax 2 Jaf_ Norb for N orb < i < N orb + N d (A15) 

corresponding to equation (A10) and equations (A13) and (A15) correspond to equation 
(All). Note that for % > N orb , —Ci is maximized at X{ = 0. The first set of Cj (for % < Norb) 
effectively keep the orbit weights positive. The second s et (for i > N nrb ) minimizes the 



components of x 2 ■ We then use the method described by iRichstone fc Tremaind (119881 ) to 
iteratively solve the equations. 

An important feature of this problem is the choice of the parameter a. Although we 
have not devised a method to specify appropriate values of this parameter in advance, if the 
quantity y 2 reflects a satisfactory accounting of the estimates of errors in the recovery of 
the LOSVDs of the target galaxy, then the choice of a and interpretation of the results is 
reasonably straightforward. A small a will generally produce a solution. Increasing a slowly 
decreases x 2 bvrt not beyond some limit. In practice, we start with a = and increase it 
slowly until the change in x 2 is l ess than 0.02 in a single iteration. 



A.3. A Test with a Black Hole 

In this case we specify a dynamical model with a gravitational potential and a distribu- 
tion function (DF) for the stars (the starlight need not, and does not, determine the mass 
in this model). The gravitational potential is 

$(r)=*4 c log(r)-^. (A16) 

r 

This potential includes a point mass and a flat circular-velocity curve (with velocity i^irc) at 
large r. We investigated Michie-like DFs of the form 

-E + J 2 /(2zut) 



f = Aexp 



a 2 



Jfn(E 1} E,E 2 ), (A17) 



where A is normalization, E and J z are energy and angular momentum, ww a (an anisotropy 
distance), N, and a (a velocity dispersion- like quantity) are model parameters, and E\ and 
E 2 are upper and lower limits on the energy. The function n is defined by 

/ „ n n X f 1 if E 1 < E < E 2 / \ \ 

"(E,E,E 2 ) = \ otlm - & " (A18) 
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The DF can be rewritten as 

/ = tiAw 2N exp(-$/a 2 ) n (v u v 2 )vf exp(-^/(2«7 2 )) exp(-t 2 /(2a 2 ))dt 2 d V(f) , (A19) 

where t 2 = v 2 ^ + v 2 , a 2 = a 2 {w 2 /{w 2 + zu 2 )) and v t = \J%E{ — This function has 
a number of interesting properties: > tends to create enhanced anisotropy and a 
flattened density distribution; w a < oo tends to depress and create a prolate density 
distribution. The cutoffs in energy avert a divergence in the luminous matter distribution 
near the black hole. 

The density is 

p = 2na 2 Auj 2N eM-®/° 2 ) / vf exp(-^/(2a 2 ))[exp(-t 2 /(2a 2 )) - exp(-t 2 /(2a 2 ))]^, 

(A20) 

where t 2 = v\ — v 2 and t\ = max(t> \ — 0). 

Note that the distribution functions with different choices of parameters can be added to 
each other (in the same potential) and similarly the densities and density-weighted moments 
are additive. A Monte Carlo realization of the model is constructed by sampling each density 
distribution and choosing the velocities randomly from the velocity distribution. Note that 
we have the additional freedom of choosing an arbitrary fraction of the velocities in the 
positive sense, and hence setting the spin of any model between a maximum amount and 
zero. The total number of random variates chosen sets the normalization. 

For the test here, we adopted a potential with f C j rc = 220 km s~ x and M. = 1.126 x 
1O 8 M . We added two DFs. The first had a = 160 km s -1 , tu a = 600 pc, N — 0, Ei — 
$(10 pc), E 2 = $(1000 pc), with equal numbers of the <fi velocities chosen in the positive 
and negative sense. The second DF had a = 120 km s -1 , w a = 200 pc, N = 2, Ex = 10 pc, 
E 2 = 1000 pc and 3/4 of the velocities were chosen positive. Each DF was assigned 
1/2 x 10 9 random points, so they had equal mass. The second DF, featuring a moderately 
spinning disky structure, dominates the model near the center of the galaxy, but declines 
rapidly. The first DF, a mildly prolate structure with no net spin dominates further out. No 
stars are found beyond a radius of 1000 pc from the center of the model. The 10 9 random 
points each are converted to projected positions and velocities (we assumed the galaxy was 
edge on), and binned into "observations" mimicing our combined HST and ground-based 
data, and also producing data mimicing data that would be obtained with an integral field 
unit. For the example below we used the IFU data with a resolution of lOpc at the center, 
and lower resolution data extending to a distance of 800 pc from the center. At each location 
the projected velocity variates were binned at a resolution of 20 km s -1 . Gaussian noise 
was added to the data, in a series of 20 realizations. These datasets were then fed into the 
modeling program, which produced (in each case) a x 2 ma P m "^circ and M.. 
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Marginalizing over each variable in turn produced the determinations of M. and t> circ 
illustrated in Figure 17. Each dataset had LOSVDs at 32 positions with 15 velocity bins, a 
product of 480. For each realization of this set, the best-fit model had a x 2 near 500. Since 
the velocity bins are uncorrelated in these datasets this is about the expected x 2 - F° r the 
presentation in Figure 17 the plots of x 2 against black hole mass and t> C i rc were arbitrarily 
shifted to a minimum of 500. The error bars judged from the individual x 2 profiles are about 
20% smaller than the error (above) determined from the spread in minima. The mean values 
of M. and t> circ from these experiments were (1.15 ± 0.03) x 1O 8 M and 222 ± 5 km s -1 , in 
excellent agreement with the target values. 

An earlier version of the code suffered from a units conversion error that led to an 
underestimate of black hole masses (and M/L). All black hole masses obtained with this 
program published in or prior to 2003 should be multiplied by 1.09. 
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Fig. I. — Ground-based NIR images of NGC 4258 identifying large-scale morphological 
characteristics of the galaxy: The bulge (PA ~ 146°), which is embedded in the inner oval 
disk (PA ~ 156°); the weak bar (PA ~ 11°); and the outer oval disk (PA ~ 150°). Top left: 
2MASS composite J+H+K image. Intensity is proportional to the square root of surface 
brightness to illustrate the bulge and the small bar. The scale is 1" pixel -1 and the image size 
is 17(1 x 17(1. Top right: As before, but intensity is proportional to surface brightness; the 
inner oval disk is now saturated and the outer oval disk becomes more prominent. Bottom 
left: K-b&nd image using SQUID on the KPNO 2.1-m telescope. Intensity is proportional 
to the square root of surface brightness to illustrate the bulge. The scale is 0"69 pixel -1 and 
the image size is 5(9 x 5(9. Bottom right: As before, but intensity is proportional to surface 
brightness to illustrate the weak bar. 




Fig. 2.— HST NICMOS imag es of NGC 4258 in the J band (F110W) (left) and in the 
K band (F222M) (right), from IChary et al.l (120001 ) . In both images, the scale and size are 
~0"038 pixel -1 and ~ 21" x 21", respectively. The images are not deconvolved. Intensity is 
proportional to the square root of surface brightness. A band of data is missing on the right 
side of the i^-band image. For the extraction of stellar mass profiles from these images, the 
dust on the SW side was corrected by replacing the SW side of the galaxy with the NE side 
flipped across the major axis. 
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Fig. 3.— HST NICMOS if- band images of th e central 5'!7 x 5'!7 of NGC 4258, adapted 
from the images provided by IChary et al.l (120001 ) . North is up, and East is to the left. Top 
left: The TinyTim model PSF shown at a nonlinear stretch (higher contrast at lower surface 
brightnesses) to emphasize the small-scale features. Top right: Same as before, but now 
scaled to the galaxy image and at the same linear stretch as in bottom row. Bottom left: 
The galaxy before PSF subtraction; the AGN is evident. Bottom right: The adopted PSF- 
subtracted galaxy image that was used to provide the central brightness profile in the K 
band. In the end, the AGN could not be subtracted well enough in K, and the J-band 
profile was used for r < 0'.'2. 



-34- 




Fig. 4. — NGC 4258 major-axis surface photometry after subtraction of prominent dust 
lanes, and correction for the AGN. The radius (r) is measured along the major axis. All 
HST profiles are deconvolved for r < 3". Top: NIR (J, J+K, K) and optical (V, R) 
surface brightness profiles. The .R-band profile is calibrated using the HST V-b&nd zeropoint 
between r ~ 5" and 7". The NIR profiles are scaled to match the optical profile at large radii 
to illustrate the absence of major color gradients in the disk. The almost power-law profile 
at r < 40" belongs to the bulge. Only data out to r w 150" are used for the dynamical 
modeling. Middle: The ellipticity as a function of radius. Bottom: The position angle of 
the major axis as a function of radius. 
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Fig. 5.— Ground-based image of NGC 4258 in the R band taken with the MDM 1.3-m 
McGraw Hill telescope. The square is centered on the nucleus, and is 8" on a side. The scale 
of the CCD was (/.'508 pixel -1 . A logarithmic stretch has been applied, and set to prevent 
"overexposure" of the bulge region. 




Fig. 6.— HST WFPC2 image of NGC 4258. The large panel shows the central 256 x 256 
pixels of the deconvolved F547M Nyquist-sampled image. It corresponds to 128 x 128 PCI 
pixels, or 5'.'84 x 5'.'84. An arbitrary logarithmic stretch has been applied. The small upper- 
right panel is the central 1'.'48 x 1'.'48 (64 x 64 subpixels) of the same image, magnified by 
a factor of 2 compared to the larger panel. Again a logarithmic stretch has been applied, 
but with greater contrast over the nuclear region. The bottom-right panel shows the same 
region divided by a model reconstructed from the surface photometry. North is up and East 
to the left. 
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Fig. 7. — Color map of the nuclear region of NGC 4258. It corresponds approximately to 
V — I, and was created from the ratio of the F547M image to the F791W image. Darker is 
redder. The map size is 8'/0 x 7'!2. The F791W image had each pixel resampled into 4 pixels 
to match the F547M image, and was then boxcar-smoothed. The central blue spike is not 
accurate because the center is saturated in the F791W data, and the F547M and F791W 
images are not perfectly aligned. 



Fig. 8. — Color composite image of NGC 4258. Red corresponds to F791W, green to F656W, 
and blue to F300W. The image has a size of 476 x 526 pixels, with a scale of 0.0455" pix" 1 . 
North is up, and East to the left. 
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Fig. 9. — Sample STIS spectra extracted from the central and two outer spatial bins on the 
approaching (SE) side of the galaxy, as well as of the template star HR 7615 which was used 
for the LOSVD deconvolution. The smooth solid lines superimposed on the galaxy spectra 
are the template stellar spectra convolved with the derived LOSVDs. The emission-line 
feature (probably Fe II) in the nuclear spectrum (dotted) was subtracted before the LOSVD 
analysis was performed. 
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Fig. 10. — As Figure [9] but for the ground-based spectra obtained with Modspec at MDM. 




Fig. 11. — Gauss-Hermite moments of the symmetrized LOSVDs as a function of position 
for the data and one model. Black filled circles are from the STIS data on the major axis. 
Red open circles are ground-based major axis data. Green triangles are ground-based minor 
axis data. The solid lines represent the analogous quantities derived from the LOSVDs of 
the base model described in 55.1. 
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Fig. 12. — NGC 4258 major-axis surface photometry: the central 10". Photometric calibra- 
tion is as in Figure HI For r < 2" there is a color gradient between the V and the NIR (J, 
J+K, and K spliced together) profiles. The V corr profile is identical to V except the inner 
three points are replaced with an arbitrary "core" profile. The NIR corr profile is created from 
NIR by replacing the points inward of r = / .'2 with points that follow the V COTT profile, shifted 
by an amount equal to the color difference between NIR and V COTT at r = 0'.'2. All profiles 
extend to r > 10" as shown in Figure HI In the NIR corr and V corr plots, only the points that 
differ from the NIR and V profiles, respectively, are shown. NIR corr and V corr are used to 
probe the effect of an unresolved AGN on the black hole mass, as discussed in §5.4[ they are 
meant to represent limiting cases, not accurate corrections to the AGN contamination. 
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Fig. 13. — Equivalent width (EW) profiles along the major axis for each component of the 
Ca II triplet (CaT) line, and for the sum of their EWs. All data points are computed from 
the STIS spectrum. The subtraction of the emission line near the center (see Figure [9]) 
may have affected the EW estimate for the 8662 A component. The 8498 A component, 
although furthest away from the emission line, has a low S/N and hence its EW profile 
may be less reliable. Furthermore, the S/N of the spectra drops rapidly with increasing 
radius. Nevertheless, the summed profile shows clearly that the spectrum within r < 0'.'2 is 
contaminated at the ~ 30% level by some source of smooth, lineless continuum light, such 
as an AGN or OB stars. The positive-r direction corresponds to the approaching (SE) side 
of the galaxy. 
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Fig. 14. — Left: Contour maps of x 2 (M., T) for the dynamical models bearing the parameters 
listed in Table HI Model names are as in Table HI The stellar M/L ratio (T) refers to the 
V band. Each dot represents a model, and dot size is proportional to the value of \ 2 f° r 
that model. The symbol "+" corresponds to the "best" model (Xmin)- The symbol "X" 
corresponds to the minimum interpolated \ 2 . The contour levels are drawn for A% 2 = 
X 2 — Xminint = 1-0 (dashed), 2.71, 4.0, and 6.63, corresponding to confidence levels of 68%, 
90%, 95%, and 99%, respectively, for one degree of freedom. The horizontal and vertical 
solid lines indicate the nominal "lex" one-dimensional uncertainties. The vertical dashed 
line indicates the maser mass. Right: Values of Xmm-mt marginalized over the T axis. The 
horizontal dashed line indicates the A% 2 = 1 level. 
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Fig. [TH — Continued 
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Fig. [TH — Continued 
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0.1 1 10 

Radius (arcsec) 

Fig. 15. — Internal moments of the baseline model in the equatorial plane. The red line 
illustrates the mean rotational velocity (v^) and the dotted line is the velocity dispersion in 
the direction a</> = J (v 2 ) — (v^) 2 . The solid line is the radial dispersion a r and the dashed 
line is erg. 
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Fig. 16. — The enclosed stellar mass M* (top row), and the sum of the black hole mass and 
the enclosed stellar mass M. + M* (bottom row), as a function of model radius. The base 
model and model P4 are constrained by the same kinematic dataset, and they share the 
same mass l^-band profile for r > 0'/04 (see Table H]). Because of the degeneracy between M. 
and M* within a small radius r ~ C.'l from the center, comparable to the spatial resolution 
of the kinematic data, the modeling algorithm can only constrain the total mass inside that 
radius. As a result, model P4 is assigned a slightly lower M. to compensate for the excess 
mass contained in the central "bump" of its mass profile. Model PI is characterized by the 
NIR corr mass profile, which differs from that of both P4 and the base model at r > r . 
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(V . /220) 2 BH mass 

Fig. 17. — x 2 profiles of orbit-based models of the Monte Carlo model described in the 
Appendix A. 3. The left panel illustrates the recovery of the circular velocity of the potential 
for 20 different realizations of its LOSVDs. The right panel illustrates the recovery of the 
black hole mass from the same set of realizations. The plots in each panel are marginalized 
over the variable in the other panel. The mean of the derived black hole masses is 1.15 x 
1O 8 M , consistent with the model mass. 
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Table 1. Journal of Imaging Observations of NGC 4258 



Instrument 


Filter/Band 


Radial Range 


Source 


HST NICMOS a 


F110W, F160W, F222M b 


r < 9" 


Charv et al. (2000) 


KPNO c 2.1-m 


K 


2" < r < 32" 


This paper 


2MASS d 


J, H, K 


r > 32" 


Jarrett et al. (2003) 


HST WFPC2 C 


F791W f 


r < 9" 


Cecil et al. (2000) 


MDM 1.3-m 


R 


r > 9" 


This paper 


HST WFPC2 


F547MS 


r < 9" 


This paper 



a Near Infrared Camera and Multi-Object Spectrometer. 

b Corresponding, approximately, to the J, H, and K bandpasses, respectively. 

c Kitt Peak National Observatory. 

d Two Micron All Sky Survey. 

e Wide Field and Planetary Camera 2. 

Corresponding, approximately, to the I bandpass. 

g Corresponding, approximately, to the V bandpass. 
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Table 2. Long-Slit Spectrograph Configurations 



Instrument a 


Slit size 


^cen 


A-range b 


Disp 1 c 


Comp. line a d 


Scale e 


+ grating 


" x " 


A 


A 


A pix -1 


A (km s" 1 ) 


" pix^ 1 


STIS + G750M 


52x0.1 


8561 


8275-8847 


0.554 


0.45 (12.4) 


0.051 


Wilbur + 831g/mm 


540x0.9 


8500 


8100-9020 


0.90 


0.9 (32) 


0.37 


Echelle + 831g/mm 


540x0.9 


8500 


8100-9020 


1.46 


0.9 (32) 


0.59 



a Wilbur and Echelle are two of the CCDs used on the Modspec instrument at the MDM 
2.4-m telescope. 

b Central wavelength and wave length range. STIS values taken from the STIS Instrument 
Handbook (ILeitherer et al.ll200ll . pp. 231, 234). For Modspec, we show the range extracted. 



c Reciprocal dispersion was measured using our own wavelength solutions. The distribution 
of dispersion solutions for a given dataset had a a ~ 0.00015 A pix . The average dispersion 
given m the Handbook for G750M is 0.56 A pix . 

instrumental line widths measured by fitting Gaussians to emission lines on comparison 
lamp exposures. This gives an estimate of the instrumental line width for extended sources. 
We use ~5 lines per exposure, and at least 5 measu rements per line . All o f the G750M ob- 
servations used an unbinned 1024x1024 pixel CCD. ILeitherer et al.l (120011 ) (p. 300) give an 
instrumental line width for point sources of a= 13.3 km s~ x for the STIS + G750M configu- 
ration. 

e The spatial scale along the slit is constant, but it varies across the slit from one grating to 
the next. It is 0'.'05597 pix~ x for G750M at 8561 A and 0"05465 pix" 1 for G750M at 6581 A 
(STIS ISR 98-23). 



Note. - - Some numbers for the CCDs Wilbur and Echelle on Modspec are calculated by 
the program "modset" by J. Thorstensen. The slit width, spectral resolution, and spatial 
resolution varied for the Modspec observations; typical measured values are shown here. 
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Table 3. Kinematic Parameter Profiles of NGC 4258 



R{") 


V 




<7 




H3 


HA 


HST STIS - Major 


■ axis profik 


j (PA = 140°) 




0.00 


20 ± 14 


159 ± 17 


0.05 ± 0.08 


-0.02 ± 0.07 


0.05 


42 ± 


6 


142 ± 


8 


-0.10 ± 0.04 


r\ r\r\ 1 f\ f\ A 

0.00 ± 0.04 


0.10 


60 ± 


4 


133 ± 


5 


-0.05 ± 0.03 


r\ r\r\ i r\ r\c\ 

-0.09 ± 0.02 


0.17 


59 ± 


4 


122 ± 


5 


-0.01 ± 0.03 


r\ r\ 'T 1 r\ r\c\ 

-0.07 ± 0.02 


0.30 


56 ± 


6 


112 ± 


7 


-0.11 ± 0.05 


C\ C\ A 1 C\ C\ A 

-0.04 ± 0.04 


0.50 


57 ± 


7 


85 ± 


6 


-0.01 ± 0.04 


-0.05 ± 0.02 


0.80 


53 ± 


8 


96 ± 10 


-0.03 ± 0.05 


-0.03 ± 0.03 


Ground-based (MDM 2.4- 


m with ModSpec) 


- Major axis profile 


0.37 


15 ± 


3 


107 ± 


6 


0.00 ± 0.04 


-0.01 ± 0.02 


0.74 


26 ± 


3 


93 ± 


6 


-0.04 ± 0.03 


n n a _i_ n no 

-(J.U4 ± (J.Uz 


1.30 


45 ± 


2 


92 ± 


5 


-0.05 ± 0.03 


-0.03 ± 0.02 


2.04 


60 ± 


2 


94 ± 


4 


-0.13 ± 0.03 


0.05 ± 0.03 


3.15 


75 ± 


3 


93 ± 


2 


-0.14 ± 0.03 


0.00 ± 0.03 


4.82 


83 ± 


3 


91 ± 


4 


-0.16 ± 0.05 


0.04 ± 0.04 


7.42 


81 ± 


5 


104 ± 


4 


-0.22 ± 0.05 


0.10 ± 0.04 


11.69 


77 ± 


7 


97 ± 


5 


-0.19 ± 0.10 


0.03 ± 0.06 


18.18 


70 ± 


8 


74 ± 


7 


-0.08 ± 0.06 


-0.05 ± 0.03 


Ground-based (MDM 2.4-: 


m with ModSpec) 


- Minor axis profile 


0.00 


-4 ± 18 


97 ± 


6 


0.02 ± 0.02 


-0.06 ± 0.01 


0.37 


1 ± 


3 


103 ± 


5 


0.01 ± 0.03 


-0.06 ± 0.02 


0.74 


3 ± 


2 


103 ± 


5 


-0.02 ± 0.04 


-0.07 ± 0.02 


1.30 


7 ± 


2 


96 ± 


4 


-0.02 ± 0.03 


-0.03 ± 0.02 


2.04 


14 ± 


2 


105 ± 


5 


0.00 ± 0.04 


0.02 ± 0.03 


3.15 


4 ± 


3 


109 ± 


4 


0.11 ± 0.03 


-0.05 ± 0.03 


4.82 


6 ± 


3 


107 ± 


4 


0.01 ± 0.04 


-0.09 ± 0.03 


7.42 


14 ± 


5 


125 ± 


6 


0.06 ± 0.06 


-0.21 ± 0.06 


11.69 


-6 ± 


5 


128 ± 


6 


0.07 ± 0.07 


-0.17 ± 0.07 


18.18 


-58 ± 11 


58 ± 14 


-0.01 ± 0.05 


-0.07 ± 0.03 



-54- 



Table 4. Dynamical Models 



Mass a Tracer 3, Kinem. b i c M. c 
Model Profile Profile Dataset (°) e d (1O 7 M ) T e 



Base 


^corr 


^corr 


S+G 


72 


0.35 


o qi +0.22 
°- o± -0.17 


3 fi +ai 
°- D -o.i 


15.11 


Dl 


^corr 


^corr 


S+G 


72 


0.45 




q 7+0.2 
°- ' -0.2 


|5.2.1| 


D2 


^corr 


^corr 


S+G 


62 


0.35 


3.621^ 


q o+O.l 


15.2.21 


Kl 


^corr 


^corr 


Sx+G 


72 


0.35 


o oc+0.21 
°- z " J -0.13 


3 fi +ai 

°-°-0.1 


|5.3| 


K2 


^corr 


^corr 


S 3 +G 


72 


0.35 


2 20 +0 - 54 

z - zu -0.31 


3 fi +ai 
°-°-o.i 


|5.3| 


K3 


Vcorr 


^corr 


G 


72 


0.35 


1 03+ 1 - 00 

1.UO_q 2g 


q y+O.l 

°- ' -0.2 


15.31 


K4 


^corr 


^corr 


S+d 


72 


0.35 


q qq + 0.21 

°- oo -0.19 


3.7^:1 


a 


PI 




^corr 


S+G 


72 


0.35 


9 40+0.09 


3.s±S:i 


I5.4| 


P2 


Vcorr 


NIR-corr 


S+G 


72 


0.35 


q qq+0.18 

o - oo -0.17 


3 fi +ai 
°-°-o.i 


15.41 


P3 


NlP^corr 


NIR, CO rr 


S+G 


72 


0.35 


2 zn+ - 21 

z -^°-0.30 


3 5 +ai 


|5.4| 


P4 






S+G 


72 


0.35 


q o c +0.22 
°- zd -0.14 


3 fi +ai 


|5.4| 




a Mass and tracer profiles 


named as 


in Fig 


;ure E2 


(cf. gOj). 







b S = All STIS kinematic data listed in Table [31 
Si = As in S but without the central kinematic datapoint. 
S3 = As in S but without the 3 central kinematic datapoints. 
G = All ground-based (Modspec) kinematic data listed in Table [3j 
Gi = As in G but without the central minor- axis datapoint. 



c Disk inclination (cf. §5.2.21) . Edge-on corresponds to i — 90°. 

d Isophote ellipticity, assumed constant throughout the galaxy (cf. §5.2.11) . 

c Black hole mass and stellar M/L ratio (T), assumed constant throughout 
the galaxy (cf. §5.ip . Uncertainty values correspond to the nominal "la" one- 
dimensional uncertainties from the contour maps (Figure [T4")) . The M/L ratio 
calibration is in the V band. 



